fourco Subroutine

      pure subroutine fourco(t,n,c,alfa,m,ress,resc,wrk1,wrk2,ier)

      !  calling sequence:
      !     call fourco(t,n,c,alfa,m,ress,resc,wrk1,wrk2,ier)
      !
      !  input parameters:
      !    t    : real array,length n, containing the knots of s(x).
      !    n    : integer, containing the total number of knots. n>=10.
      !    c    : real array,length n, containing the b-spline coefficients.
      !    alfa : real array,length m, containing the parameters alfa(i).
      !    m    : integer, specifying the number of integrals to be computed.
      !    wrk1 : real array,length n. used as working space
      !    wrk2 : real array,length n. used as working space
      !
      !  output parameters:
      !    ress : real array,length m, containing the integrals ress(i).
      !    resc : real array,length m, containing the integrals resc(i).
      !    ier  : error flag:
      !      ier=0 : normal return.
      !      ier=10: invalid input data (see restrictions).
      !
      !  restrictions:
      !    n >= 10
      !    t(4) < t(5) < ... < t(n-4) < t(n-3).
      !    t(1) <= t(2) <= t(3) <= t(4).
      !    t(n-3) <= t(n-2) <= t(n-1) <= t(n).
      !
      !  other subroutines required: fpbfou,fpcsin
      !
      !  references :
      !    dierckx p. : calculation of fourier coefficients of discrete functions using cubic splines.
      !                 j. computational and applied mathematics 3 (1977) 207-209.
      !    dierckx p. : curve and surface fitting with splines, monographs on
      !                 numerical analysis, oxford university press, 1993.
      !
      !  author :
      !    p.dierckx
      !    dept. computer science, k.u.leuven
      !    celestijnenlaan 200a, b-3001 heverlee, belgium.
      !    e-mail : Paul.Dierckx@cs.kuleuven.ac.be
      !
      !  ..scalar arguments..
      integer, intent(in) :: n,m
      integer, intent(out) :: ier
      !  ..array arguments..
      real(RKIND), intent(in) :: t(n),c(n),alfa(m)
      real(RKIND), intent(inout) :: wrk1(n),wrk2(n)
      real(RKIND), intent(out) :: ress(m),resc(m)
      !  ..local scalars..
      integer :: i,n4
      !  ..
      n4 = n-4
      !  before starting computations a data check is made. in the input data
      !  are invalid, control is immediately repassed to the calling program.
      ier = FITPACK_INPUT_ERROR

      ! Not enough points
      if (n<10) return

      ! Ends of the support: knots must be monotonic
      if (any(t(1:3)>t(2:4)))       return
      if (any(t(n-2:n)<t(n-3:n-1))) return

      ! Interior: knots must be strictly monotonic
      if (any(t(4:n4)>=t(5:n-3)))   return

      ier = FITPACK_OK
      !  main loop for the different alfa(i).
      alphas: do i=1,m

         !  calculate the integrals
         !    wrk1(j) = integral(nj,4(x)*sin(alfa*x))    and
         !    wrk2(j) = integral(nj,4(x)*cos(alfa*x)),  j=1,2,...,n-4,
         !  where nj,4(x) denotes the normalised cubic b-spline defined on the knots t(j),t(j+1),...,t(j+4).
         call fpbfou(t,n,alfa(i),wrk1,wrk2)

         ! calculate the integrals ress(i) and resc(i).
         ress(i) = dot_product(c(1:n4),wrk1(1:n4))
         resc(i) = dot_product(c(1:n4),wrk2(1:n4))
      end do alphas

      end subroutine fourco